####################################################################################
# Elections Improve Support for State Trial Court Judges in the United States ######
# 2021 Survey Analyses #############################################################
# Andrew R. Stone and Michael P. Olson #############################################
####################################################################################

####################################################################################
# Set working directory, load packages, functions ##################################
####################################################################################
  # Set working directory
  # Load packages
  library(cjoint); library(ggcorrplot); library(ggplot2); library(psych); library(questionr); library(sandwich); library(stargazer); library(tidyverse); library(vtable)

  # Prediction and robust SE function
  predict.rob <- function(x,clcov,newdata){
    if(missing(newdata)){ newdata <- x$model }
    tt <- terms(x)
    Terms <- delete.response(tt)
    m.mat <- model.matrix(Terms,data=newdata)
    m.coef <- x$coef
    fit <- as.vector(m.mat %*% x$coef)
    se.fit <- sqrt(diag(m.mat%*%clcov%*%t(m.mat)))
    return(list(fit=fit,se.fit=se.fit))}
  
  # Ceiling function for p values
  p_ceiling <- function(x){return(ceiling(100*x)/100)}

#############
# Load data #
#############
  
  load("weidenbaum_survey_2021.RData")

########################################################################################
# Descriptive Statistics ###############################################################
########################################################################################
#############
# Table A.2 #
#############
# Summary statistics for outcome variables
  
  # Support for judge
  # N
  cat(nrow(survey[!is.na(survey$judge_support),]), sep = '\n', file = "./figures_tables/n_support.tex")
  # Values
  cat(paste(as.character(format(round(prop.table(table(survey$judge_support)),2),nsmall=2)),collapse=" & ",sep=" & "), 
      sep = '\n', file = "./figures_tables/row_support.tex")
  
  # Judge fairness
  # N
  cat(nrow(survey[!is.na(survey$fairness_accused_crime),]), sep = '\n', file = "./figures_tables/n_fairness.tex")
  # Values
  cat(paste(as.character(format(round(prop.table(table(survey$fairness_accused_crime)),2),nsmall=2)),collapse=" & ",sep=" & "), 
      sep = '\n', file = "./figures_tables/row_fairness.tex")
  
  # Legitimacy
  # N (additive)
  cat(nrow(survey[!is.na(survey$legitimacy_additive),]), sep = '\n', file = "./figures_tables/n_legitimacy.tex")
  # Values (individual questions)
  # These are flipped because they are reverse coded to make higher values more legitimate
  cat(paste(as.character(format(round(prop.table(table(6-survey$legitimacy_1)),2),nsmall=2)),collapse=" & ",sep=" & "), 
      sep = '\n', file = "./figures_tables/row_doaway.tex")
  cat(paste(as.character(format(round(prop.table(table(6-survey$legitimacy_2)),2),nsmall=2)),collapse=" & ",sep=" & "), 
      sep = '\n', file = "./figures_tables/row_remove.tex")
  cat(paste(as.character(format(round(prop.table(table(6-survey$legitimacy_3)),2),nsmall=2)),collapse=" & ",sep=" & "), 
      sep = '\n', file = "./figures_tables/row_independent.tex")
  cat(paste(as.character(format(round(prop.table(table(6-survey$legitimacy_4)),2),nsmall=2)),collapse=" & ",sep=" & "), 
      sep = '\n', file = "./figures_tables/row_mixed.tex")

#############
# Table A.1 #
#############
# Summary statistics for respondent characteristics

  # Experience with elections (for all states, if elections, trial courts are elected)
  survey$Experience <- recode(survey$elected_any_court,"1"="Any Court Elected in State","0"="No Courts Elected in State")
  # Gender/sex
  survey$GENDER <- recode(survey$GENDER,"1"="Male","2"="Female")
  survey$Sex <- factor(survey$GENDER,levels=c("Female","Male"))
  # Race/ethnicity (3 is other, 5 is multiple, we combine)
  survey$RACETHNICITY <- recode(survey$RACETHNICITY,"1"="White","2"="Black",
                                "3"="Other/Multiple","4"="Hispanic","5"="Other/Multiple",
                                "6"="Asian")
  survey$`Race or Ethnicity` <- factor(survey$RACETHNICITY)
  # Partisanship
  survey$Party <- factor(survey$respondent_pid)
  # Education level
  survey$EDUC5 <- recode(survey$EDUC5,"1"="Less than High School","2"="High School or Equivalent",
                                "3"="Vocational School or Some College","4"="Bachelor's Degree","5"="Post-Graduate Study")
  survey$Education <- factor(survey$EDUC5,
                             levels=c("Less than High School",
                                      "High School or Equivalent",
                                      "Vocational School or Some College",
                                      "Bachelor's Degree",
                                      "Post-Graduate Study"))
  # Income level
  survey$INCOME4 <- recode(survey$INCOME4,"1"="Less than $30,000","2"="$30,000 to $60,000",
                         "3"="$60,000 to $100,000","4"="$100,000 or More")
  survey$Income <- factor(survey$INCOME4,
                          levels=c("Less than $30,000",
                                   "$30,000 to $60,000",
                                   "$60,000 to $100,000",
                                   "$100,000 or More"))
  # Ideology
  survey$IDEO <- recode(survey$IDEO,"1"="Very Liberal","2"="Somewhat Liberal",
                           "3"="Moderate","4"="Somewhat Conservative","5"="Very Conservative")
  survey$Ideology <- factor(survey$IDEO,
                            levels=c("Very Liberal",
                                     "Somewhat Liberal",
                                     "Moderate",
                                     "Somewhat Conservative",
                                     "Very Conservative"))
  
  # Summary statistics table with these variables
  sumtab <- sumtable(data=survey,
           vars=c("Sex","Race or Ethnicity","Party",
                  "Education","Income","Ideology","Experience"),
           out="return",
           title="Descriptive Statistics: July 2021 NORC Survey")
  # Feeding into stargazer
  sumtab <- stargazer(sumtab,summary=F,
            title="Descriptive Statistics: July 2021 NORC Survey",
            label="norc-demographics",rownames = F)
  # Formatting for LaTeX
  sumtab <- gsub("ccc","lcc",sumtab,fixed=T)
  sumtab <- gsub("Sex","\\\\textbf{Sex}",sumtab)
  sumtab <- gsub("Race or Ethnicity","\\\\textbf{Race or Ethnicity}",sumtab)
  sumtab <- gsub("Party","\\\\textbf{Party}",sumtab)
  sumtab <- gsub("Education","\\\\textbf{Education}",sumtab)
  sumtab <- gsub("Income","\\\\textbf{Income}",sumtab)
  sumtab <- gsub("Ideology","\\\\textbf{Ideology}",sumtab)
  sumtab <- gsub("Experience","\\\\textbf{Experience}",sumtab)
  # Saving the table
  cat(sumtab, sep = '\n', file = "./figures_tables/studyonesumtab.tex")

############################################
# In-text discussion of legitimacy measure #
############################################
# Cronbach's alpha for legitimacy
  # Estimating with the four legitimacy questions
  alpha <- alpha(data.frame(survey$legitimacy_1,survey$legitimacy_2,survey$legitimacy_3,survey$legitimacy_4))
  # Extracting the estimate and saving
  cat(as.character(format(round(alpha$total[1],2),nsmall=2)), sep = '\n', file = "./figures_tables/legitimacy_alpha.tex")

##############
# Figure A.1 #
##############
# Correlation matrix of outcome variables
  # Outcome variables to make the matrix out of
  dat_bit <- survey[,c("judge_support","fairness_accused_crime",
                       "legitimacy_1","legitimacy_2","legitimacy_3","legitimacy_4")]
  # Keeping only respondents who answered all outcome questions
  dat_bit <- dat_bit[complete.cases(dat_bit),]
  # Renaming columns for clarity
  colnames(dat_bit) <- c("Support",
                         "Fairness",
                         "Do Away with Court",
                         "Remove Judges",
                         "Courts too Independent",
                         "Mixed Up in Politics")
  # Creating the correlation matrix
  ggcorrplot(cor(dat_bit),
             outline.color = "white",
             type="lower",show.legend = F,
             lab=T,
             colors=c("indianred","white","grey20"))
  # Saving the correlation matrix
  
  ggsave("./figures_tables/cormat.eps",width=8,height=6,device=cairo_ps)
  # Saving the number of respondents who answered all outcome questions
  cat(as.character(format(nrow(dat_bit),big.mark=",")), sep = '\n', file = "./figures_tables/obs_cors.tex")
  
##########################################################################
# Treatment Balance ######################################################
##########################################################################  
#############
# Table B.1 #
#############  
  # Using OLS, regressing respondent characteristics on profile attributes
  # Limiting to respondents who got treatment assignments
  survey.name <- survey[!is.na(survey$judge_name),]
  
  black_outcome <- lm(black ~ judge_name + judge_race 
                      + election_treatment_binary + judge_margin 
                      + judge_tenure + judge_party 
                      + judge_sentencing, 
                      data=survey.name)
  
  hispanic_outcome <- lm(hispanic ~ judge_name + judge_race 
                         + election_treatment_binary + judge_margin 
                         + judge_tenure + judge_party 
                         + judge_sentencing, 
                         data=survey.name)
  
  college_outcome <- lm(college ~ judge_name + judge_race 
                        + election_treatment_binary + judge_margin 
                        + judge_tenure + judge_party 
                        + judge_sentencing, 
                        data=survey.name)
  
  married_outcome <- lm(married ~ judge_name + judge_race 
                        + election_treatment_binary + judge_margin 
                        + judge_tenure + judge_party 
                        + judge_sentencing, 
                        data=survey.name)
  
  employed_outcome <- lm(employed ~ judge_name + judge_race 
                         + election_treatment_binary + judge_margin 
                         + judge_tenure + judge_party 
                         + judge_sentencing, 
                         data=survey.name)
  
  income_outcome <- lm(INCOME ~ judge_name + judge_race 
                       + election_treatment_binary + judge_margin 
                       + judge_tenure + judge_party 
                       + judge_sentencing, 
                       data=survey.name)
  
  internet_outcome <- lm(INTERNET ~ judge_name + judge_race 
                         + election_treatment_binary + judge_margin 
                         + judge_tenure + judge_party 
                         + judge_sentencing, 
                         data=survey.name)
  
  democrat_outcome <- lm(democrat ~ judge_name + judge_race 
                         + election_treatment_binary + judge_margin 
                         + judge_tenure + judge_party 
                         + judge_sentencing, 
                         data=survey.name)
  
  republican_outcome <- lm(republican ~ judge_name + judge_race 
                           + election_treatment_binary + judge_margin 
                           + judge_tenure + judge_party 
                           + judge_sentencing, 
                           data=survey.name)
  # List of models
  balance_models_list <- list(black_outcome,hispanic_outcome,college_outcome,married_outcome,
                              employed_outcome,income_outcome,internet_outcome,democrat_outcome,
                              republican_outcome)
  # Stargazer table of output
  balance_sg <- stargazer(balance_models_list,
                          se=list(sqrt(diag(vcovHC(black_outcome,type="HC2"))), # HC2 matches the output of amce (the cjoint function)
                                  sqrt(diag(vcovHC(hispanic_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(college_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(married_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(employed_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(income_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(internet_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(democrat_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(republican_outcome,type="HC2")))
                          ),
                          dep.var.labels = c("Black","Hispanic","College","Married","Employed",
                                             "Income","Internet","Democrat","Republican"),
                          covariate.labels = c("Woman","White","Hispanic","Appointed",
                                               "60\\% Margin","80\\% Margin",
                                               "5 Year Tenure","15 Year Tenure",
                                               "Democrat","Party Not Shown","6 Year Avg. Sentence",
                                               "9 Year Avg. Sentence","Constant"),
                          keep.stat=c("n","adj.rsq"),
                          notes.append = FALSE,
                          star.char="*",star.cutoffs=0.05,notes.label = "",
                          font.size="scriptsize",column.sep.width="0cm",no.space=T,
                          notes="\\parbox[t]{0.95\\textwidth}{\\scriptsize \\textit{Note}: Estimates are from OLS regressions with robust standard errors
                        in parentheses. The unit of observation is the respondent. All outcomes are binary except 
                        income, which is measured on an 18-point scale and is treated as a continuous variable. $^{*}p<0.05$.}",
                          title="Balance Test Results",
                          label="balance"
  )
  # Saving the table
  cat(balance_sg, sep = '\n', file = "./figures_tables/balance.tex")

####################################################################################
# Conjoint analyses ################################################################
####################################################################################
#######################  
# Set up the conjoint #
#######################
  # Names of attributes to plot
  attribute.names.full <- c("Selection Method",
                            "Gender","Race/Ethnicity",
                            "Selection Margin","Tenure",
                            "Partisanship","Sentencing Behavior")
  
  # Names of varying attributes to plot
  levels.test.full<-list()
  levels.test.full[["a.election_treatment"]]<-c("Appointed","Partisan election","Nonpartisan election") # Baseline category is set below
  levels.test.full[["b.judge_name"]]<-c("Man","Woman")
  levels.test.full[["c.judge_race"]]<-c("Black","White","Hispanic")
  levels.test.full[["d.judge_margin"]]<-c("52","60","80")
  levels.test.full[["e.judge_tenure"]]<-c("1","5","15")
  levels.test.full[["f.judge_party"]]<-c("Republican","Democrat","Not Listed")
  levels.test.full[["g.judge_sentencing"]]<-c("Average","Below","Above") # Baseline category is set below
  
#########################
# Support for the judge #
#########################
############
# Figure 1 #
############
  # Variables to go in the model, ensuring they are alphabetically ordered in way we wish to present them
  survey$a.election_treatment <- survey$election_treatment
  survey$b.judge_name <- survey$judge_name
  survey$c.judge_race <- survey$judge_race
  survey$d.judge_margin <- survey$judge_margin
  survey$e.judge_tenure <- survey$judge_tenure
  survey$f.judge_party <- survey$judge_party
  survey$g.judge_sentencing <- survey$judge_sentencing
  
  # Adjusting baselines so sentencing baseline is 6 years (the average) and the election treatment is appointed
  baselines_list <- list()
  baselines_list$g.judge_sentencing <- "2"
  baselines_list$a.election_treatment <- "3"
  
  # Subset dataset to respondents who answered the support question 
  survey.support <- survey[!(is.na(survey$judge_support)),]
  survey.support$judge_support <- as.numeric(as.character(survey.support$judge_support))
  
  # Estimate the basic AMCEs, all election types
  results.support <- amce(judge_support ~ a.election_treatment + b.judge_name  + c.judge_race + 
                                          d.judge_margin + e.judge_tenure + f.judge_party + g.judge_sentencing, 
                                                      data=survey.support, baselines=baselines_list, cluster=F)

  # Create figure
  plot(results.support, attribute.names=attribute.names.full,
                      label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1.25),
                      color="black", ylim=c(-1,1))
  ggsave("./figures_tables/cjoint-support.eps", width=7, height=5,device=cairo_ps)

######################################################
# In-text discussion of Figure 1 and related results #
######################################################
  
  # Extracting coefficients and p-values to discuss in the text
  # OLS version replicating results from above (including HC2 robust standard errors)
  results.support.ols <- lm(judge_support ~ judge_name + judge_race + I(election_treatment==1) + I(election_treatment==2) + judge_margin + judge_tenure 
                          + judge_party + judge_sentencing, 
                          data=survey.support)
  results.support.ols <- coeftest(results.support.ols,vcovHC(results.support.ols,type="HC2"))
  # Coefficient and p-value on the partisan election treatment
  cat(round(results.support.ols["I(election_treatment == 1)TRUE","Estimate"],2), sep = '\n', file = "./figures_tables/partisan_coef.tex")
  cat(p_ceiling(results.support.ols["I(election_treatment == 1)TRUE","Pr(>|t|)"]), sep = '\n', file = "./figures_tables/partisan_p.tex")
  # Coefficient and p-value on the nonpartisan election treatment
  cat(round(results.support.ols["I(election_treatment == 2)TRUE","Estimate"],2), sep = '\n', file = "./figures_tables/nonpartisan_coef.tex")
  cat(p_ceiling(results.support.ols["I(election_treatment == 2)TRUE","Pr(>|t|)"]), sep = '\n', file = "./figures_tables/nonpartisan_p.tex")
  # Sentencing behavior (below) coefficient and p-value
  round(summary(results.support)$amce$Estimate[12],2)
  p_ceiling(summary(results.support)$amce$"Pr(>|z|)"[12])
  # Sentencing behavior (above) coefficient and p-value
  round(summary(results.support)$amce$Estimate[13],2)
  p_ceiling(summary(results.support)$amce$"Pr(>|z|)"[13])
  
  # Analysis of binary election treatment
  # OLS, mirroring cjoint analysis (including HC2 robust standard errors)
  results.support.binary <- lm(judge_support ~ judge_name + judge_race + I(election_treatment_binary==1) + judge_margin + judge_tenure 
                               + judge_party + judge_sentencing, 
                               data=survey.support)
  results.support.binary <- coeftest(results.support.binary,vcovHC(results.support.binary,type="HC2"))
  # Coefficient and p-value on the binary election treatment
  cat(round(results.support.binary["I(election_treatment_binary == 1)TRUE","Estimate"],2), sep = '\n', file = "./figures_tables/binary_coef.tex")
  cat(p_ceiling(results.support.binary["I(election_treatment_binary == 1)TRUE","Pr(>|t|)"]), sep = '\n', file = "./figures_tables/binary_p.tex")
  
  # Comparison of average magnitude of sentencing effect and binary election treatment 
  mean(c(round(summary(results.support)$amce$Estimate[12],2), round(summary(results.support)$amce$Estimate[13],2)))/round(results.support.binary["I(election_treatment_binary == 1)TRUE","Estimate"],2)
  
########################################################
# In-text discussion of shared race and gender results #
########################################################
  # Coding shared race
  survey.support$shared_race <- 0
  survey.support[survey.support$judge_race==1 & survey.support$`Race or Ethnicity`=="Black","shared_race"] <- 1
  survey.support[survey.support$judge_race==3 & survey.support$`Race or Ethnicity`=="Hispanic","shared_race"] <- 1
  survey.support[survey.support$judge_race==2 & survey.support$`Race or Ethnicity`=="White","shared_race"] <- 1
  # T-test of difference in support for shared race and not shared race
  shared_race <- t.test(survey.support$judge_support[survey.support$shared_race==1],
                        survey.support$judge_support[survey.support$shared_race==0])
  # p-value on effect of shared race
  cat(as.character(format(nsmall=2,p_ceiling(shared_race$p.value))), sep = '\n', file = "./figures_tables/shared_race_p.tex")
  
  # Coding shared gender
  survey.support$shared_gender <- 0
  survey.support[survey.support$Sex=="Male" & survey.support$judge_name==1,"shared_gender"] <- 1
  survey.support[survey.support$Sex=="Female" & survey.support$judge_name==2,"shared_gender"] <- 1
  # T-test of difference in support for shared gender and not shared gender
  shared_gender <- t.test(survey.support$judge_support[survey.support$shared_gender==1],
                          survey.support$judge_support[survey.support$shared_gender==0])
  # p-value on effect of shared gender
  cat(as.character(format(nsmall=2,p_ceiling(shared_gender$p.value))), sep = '\n', file = "./figures_tables/shared_gender_p.tex")
  
##############################################
# Support for the judge as function of party #
##############################################
############
# Figure 2 #
############
  # Estimate the basic AMCEs, varying by respondent party
  results.support.by.party <- amce(judge_support ~ 
                                     a.election_treatment*Party + b.judge_name*Party  + c.judge_race*Party + 
                                     d.judge_margin*Party + e.judge_tenure*Party + f.judge_party*Party + g.judge_sentencing*Party, 
                          data=survey.support, respondent.varying="Party", baselines = baselines_list, cluster=F,
                          na.ignore = T)
  # Create figure
  plot(results.support.by.party, plot.display="interaction", attribute.names = attribute.names.full,
       label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1.25),
       color="black")

  ggsave("./figures_tables/cjoint-support-by-party.eps", width=7, height=5,device=cairo_ps)
  
  # Size of subgroups for each analysis
  table(survey.support$Party)
  
###############################################
# In-text discussion of results from Figure 2 #
###############################################
  
  # Democratic respondent support for Democratic judges as opposed to Republican
  round(summary(results.support.by.party)$Party1amce$Estimate[10],2)
  p_ceiling(summary(results.support.by.party)$Party1amce$"Pr(>|z|)"[10])
  # Approximate increase in support (scale is from 1-5, so dividing by 4)
  (round(summary(results.support.by.party)$Party1amce$Estimate[10],2)/4)*100
  
  # Republican respondent support for Republican judges as opposed to Democratic
  round(summary(results.support.by.party)$Party3amce$Estimate[10],2)*-1
  p_ceiling(summary(results.support.by.party)$Party3amce$"Pr(>|z|)"[10])
  # Approximate increase in support (scale is from 1-5, so dividing by 4)
  (round(summary(results.support.by.party)$Party3amce$Estimate[10],2)/4)*-1*100
  
  # Comparison to magnitude of combined selection method effect from above model
  mean(c((round(summary(results.support.by.party)$Party3amce$Estimate[10],2))*-1, round(summary(results.support.by.party)$Party1amce$Estimate[10],2)))/round(results.support.binary["I(election_treatment_binary == 1)TRUE","Estimate"],2)
  
#######################################################################  
# Comparison of treatment effect of elections for each partisan group #
#######################################################################  
  # Calculating difference on election treatment estimate for I, R compared to D
  # I compared to D
  # Partisan (p value from z score, which is estimate/SE)
  2*pnorm(q=results.support.by.party$cond.estimates$`aelectiontreatment:Party`[1,1]/results.support.by.party$cond.estimates$`aelectiontreatment:Party`[2,1], lower.tail = F)
  # Nonpartisan
  2*pnorm(q=results.support.by.party$cond.estimates$`aelectiontreatment:Party`[1,2]/results.support.by.party$cond.estimates$`aelectiontreatment:Party`[2,2], lower.tail = F)

  # R compared to D
  # Partisan
  2*pnorm(q=results.support.by.party$cond.estimates$`aelectiontreatment:Party`[1,3]/results.support.by.party$cond.estimates$`aelectiontreatment:Party`[2,3], lower.tail = F)
  # Nonpartisan
  2*pnorm(q=results.support.by.party$cond.estimates$`aelectiontreatment:Party`[1,4]/results.support.by.party$cond.estimates$`aelectiontreatment:Party`[2,4], lower.tail = F)

  # Saving for in-text discussion of p-value on difference of partisan treatment for I and D 
  i_vs_d_part <- 2*pnorm(q=results.support.by.party$cond.estimates$`aelectiontreatment:Party`[1,1]/results.support.by.party$cond.estimates$`aelectiontreatment:Party`[2,1], lower.tail = F)
  cat(as.character(format(p_ceiling(i_vs_d_part))), sep = '\n', file = "./figures_tables/p_i_vs_d_partisan.tex")
  
  # Saving for in-text discussion of treatment effect and p-value on difference of nonpartisan treatment for R and D 
  r_vs_d_nonpart <- 2*pnorm(q=results.support.by.party$cond.estimates$`aelectiontreatment:Party`[1,4]/results.support.by.party$cond.estimates$`aelectiontreatment:Party`[2,4], lower.tail = F)
  cat(as.character(format(round(results.support.by.party$cond.estimates$`aelectiontreatment:Party`[1,4],2),nsmall=2)), sep = '\n', file = "./figures_tables/coef_r_vs_d_nonpartisan.tex")
  cat(as.character(format(p_ceiling(r_vs_d_nonpart),nsmall=2)), sep = '\n', file = "./figures_tables/p_r_vs_d_nonpartisan.tex")

  # Version of analysis with I as baseline (to readily compare I and R)
  survey.support$Party_I_baseline <- relevel(survey.support$Party, ref="Independent")
  results.support.by.party.i.baseline <- amce(judge_support ~ 
                                   a.election_treatment*Party_I_baseline + b.judge_name*Party_I_baseline  + c.judge_race*Party_I_baseline + 
                                   d.judge_margin*Party_I_baseline + e.judge_tenure*Party_I_baseline + f.judge_party*Party_I_baseline + g.judge_sentencing*Party_I_baseline, 
                                 data=survey.support, respondent.varying="Party_I_baseline", baselines = baselines_list, cluster=F,
                                 na.ignore = T)

  # I compared to R (I compared to D we get from the model above)
  # Partisan (p value from z score, which is estimate/SE)
  2*pnorm(q=(results.support.by.party.i.baseline$cond.estimates$`aelectiontreatment:PartyIbaseline`[1,3]/results.support.by.party.i.baseline$cond.estimates$`aelectiontreatment:PartyIbaseline`[2,3])*-1, lower.tail = F)
  # Nonpartisan (p value from z score, which is estimate/SE)
  2*pnorm(q=(results.support.by.party.i.baseline$cond.estimates$`aelectiontreatment:PartyIbaseline`[1,4]/results.support.by.party.i.baseline$cond.estimates$`aelectiontreatment:PartyIbaseline`[2,4]), lower.tail = F)

#####################################################################
# Evaluation of judge fairness in hearing case of someone they knew #
#####################################################################
##############
# Figure B.1 #
##############
  # Subset dataset to respondents who answered the judge fairness question
  survey.judge.fairness <- survey[!(is.na(survey$fairness_accused_crime)),]
  survey.judge.fairness$fairness_accused_crime <- as.numeric(survey.judge.fairness$fairness_accused_crime)

  # Estimate the basic AMCEs, fairness
  results.fairness <- amce(fairness_accused_crime ~ a.election_treatment + b.judge_name  + c.judge_race + 
                             d.judge_margin + e.judge_tenure + f.judge_party + g.judge_sentencing, 
                           data=survey.judge.fairness, cluster = F , baselines = baselines_list)

# Create figure
  
  plot(results.fairness, attribute.names = attribute.names.full,
                          label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1.25),
                          color="black")
  ggsave("./figures_tables/cjoint-fairness.eps", width=7, height=5,device=cairo_ps)

##############
# Figure B.2 #
##############
  # Estimate the AMCEs by party, fairness
  results.fairness.by.party <- amce(fairness_accused_crime ~ a.election_treatment*Party + b.judge_name*Party  + c.judge_race*Party + 
                                      d.judge_margin*Party + e.judge_tenure*Party + f.judge_party*Party + g.judge_sentencing*Party, na.ignore = T,
                                    data=survey.judge.fairness, respondent.varying="Party", cluster = F , baselines = baselines_list)

# Create figure
  
  plot(results.fairness.by.party, plot.display="interaction", attribute.names = attribute.names.full,
                                   label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1),
                                   color="black")

  ggsave("./figures_tables/cjoint-fairness-by-party.eps", width=7, height=5,device=cairo_ps)

#########################
# Legitimacy evaluation #
#########################
############
# Figure 3 #
############
  # Subset dataset to respondents who answered the legitimacy questions
  survey.legitimacy <- survey[!(is.na(survey$legitimacy_additive)),]

# Estimate the basic AMCEs, legitimacy
  results.legitimacy <- amce(legitimacy_additive ~ a.election_treatment + b.judge_name  + c.judge_race + 
                               d.judge_margin + e.judge_tenure + f.judge_party + g.judge_sentencing,
                             data=survey.legitimacy,cluster = F, baselines = baselines_list)
  
# Create figure

  figure.legitimacy <- plot(results.legitimacy, attribute.names = attribute.names.full,
                            label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1.25),
                            color="black")
  dev.off()
  ggsave("./figures_tables/cjoint-legitimacy.eps", width=7, height=5,device=cairo_ps)
  
##############
# Figure B.3 #
##############
  # Legitimacy analysis, removing the state high court question
  # Subset dataset to respondents who answered the legitimacy questions (removing state high court question)
  survey.legitimacy.nosc <- survey[!(is.na(survey$legitimacy_nosc)),]
  
  # Estimate the basic AMCEs, legitimacy, no state high court question
  results.legitimacy.nosc <- amce(legitimacy_nosc ~ a.election_treatment + b.judge_name  + c.judge_race + 
                               d.judge_margin + e.judge_tenure + f.judge_party + g.judge_sentencing,
                             data=survey.legitimacy.nosc,cluster = F, baselines = baselines_list)
  
  # Create figure
  figure.legitimacy.nosc <- plot(results.legitimacy.nosc, attribute.names = attribute.names.full,
                            label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1.25),
                            color="black")
  ggsave("./figures_tables/cjoint-legitimacy-nosc.eps", width=7, height=5,device=cairo_ps)

###################################
# Legitimacy as function of party #
###################################
##############
# Figure B.4 #
##############
  # Estimate the AMCEs by party, legitimacy
  results.legitimacy.by.party <- amce(legitimacy_additive ~ a.election_treatment*Party + b.judge_name*Party  + c.judge_race*Party + 
                                        d.judge_margin*Party + e.judge_tenure*Party + f.judge_party*Party + g.judge_sentencing*Party, na.ignore = T,
                                      data=survey.legitimacy, respondent.varying="Party",cluster = F, baselines = baselines_list)

  # Create figure
  figure.legitimacy.by.party <- plot(results.legitimacy.by.party, facet.names = "Party", plot.display="interaction", attribute.names = attribute.names.full,
                                     label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1,1),
                                     color="black")
  ggsave("./figures_tables/cjoint-legitimacy-by-party.eps", width=7, height=5,device=cairo_ps)

##################################################
# Support for the judge as function of knowledge #
##################################################
##############
# Figure B.5 #
##############
  # Estimate the basic AMCEs by knowledge
  
  survey.support$`Knowledge Level` <- factor(survey.support$`Knowledge Level`,levels=c("Not Knowledgeable","Knowledgeable"))
  
  results.support.by.smart <- amce(judge_support ~ a.election_treatment*`Knowledge Level` + b.judge_name*`Knowledge Level` + c.judge_race*`Knowledge Level` + 
                                     d.judge_margin*`Knowledge Level` + e.judge_tenure*`Knowledge Level` + f.judge_party*`Knowledge Level` + g.judge_sentencing*`Knowledge Level`,na.ignore = T,
                                   data=survey.support, respondent.varying="Knowledge Level", baselines = baselines_list, cluster=F)
  
  # Create figure
  plot(results.support.by.smart, plot.display="interaction", attribute.names = attribute.names.full,
                                  label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1.25),
                                  color="black")

  ggsave("./figures_tables/cjoint-support-by-knowledge.eps", width=7, height=5,device=cairo_ps)

##########################################
# Interactive effects, support for judge #
##########################################
#######################
# Figure E.1, Panel a #
#######################
  # Subset to respondents who are partisans and got a partisan judge treatment
  survey.support.partisans <- survey.support[which(survey.support$respondent_pid %in% c("Democrat","Republican") & survey.support$judge_party!=3),]
  # Copartisan variable
  survey.support.partisans$judge_copartisan <- ifelse((survey.support.partisans$judge_party == 2 & survey.support.partisans$respondent_pid == "Democrat") |
                                                        (survey.support.partisans$judge_party == 1 & survey.support.partisans$respondent_pid == "Republican"), 1, 0)

  # OLS, support on indicator for the binary election treatment interacted with copartisanship
  copartisan_reg <- lm(judge_support ~ judge_copartisan*election_treatment_binary, data=survey.support.partisans)
  
  # Predicted values
  copartisan_data <- data.frame(election_treatment_binary=as.factor(c(1,2,1,2)), judge_copartisan=c(0,0,1,1))
  copartisan_fitted <- predict.rob(x=copartisan_reg,newdata=copartisan_data,clcov=vcovHC(copartisan_reg,type="HC0"))

  # Creating the plot
  coef <- c(copartisan_fitted$fit[1],copartisan_fitted$fit[3],copartisan_fitted$fit[2],copartisan_fitted$fit[4])
  se <- c(copartisan_fitted$se.fit[1],copartisan_fitted$se.fit[3],copartisan_fitted$se.fit[2],copartisan_fitted$se.fit[4])
  lb <- coef-1.96*se
  ub <- coef+1.96*se
  var <- c("Elected","Elected","Appointed","Appointed")
  party <- c("Outpartisan","Copartisan","Outpartisan","Copartisan")

  results <- data.frame(coef,se,lb,ub,var,party)
  results$var <- reorder(results$var,rep(1:2,each=2))
  
  ggplot(results, aes(x=coef,y=var,color=party,group=party)) +
    geom_point(size=4) + 
    geom_segment(aes(x=lb,xend=ub,y=var,yend=var),size=1.25)+
    scale_color_grey(end=0.6)+
    theme_bw()+
    xlab("Predicted Support")+
    labs(colour="")+
    theme(legend.position="bottom",
          axis.title.y = element_blank())
  
  # Saving the plot
  ggsave("./figures_tables/elections-moderate-differences.eps",width=6,height=4,device=cairo_ps)
  